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I.  INTRODUCTION 


In  recent  years,  a  considerable  research  effort  has  been  focused  on  the 
development  of  modern  predictive  capabilities  for  determining  the  aerodynamics 
of  projectiles.  The  time-dependent  Navier-Stokes  computational  technique  has 
been  used1”2  to  compute  the  flow  over  projectiles  at  transonic  speeds.  For 
supersonic  flows,  space-marching  parabolized3  Navier-Stokes  computational 
technique  can  be  effectively  used.  However,  this  technique  fails  for  flows 
containing  longitudinal  flow  separation.  In  such  cases,  which  are  frequently 
encountered  in  projectile  aerodynamic  simulations,  the  time-dependent  Navier- 
Stokes  technique  needs  to  be  used.4”5 

The  time-dependent  Navier-Stokes  equations  can  be  solved  in  a  generalized 
body-fitted  coordinate  system.  Many  actual  projectile  configurations  contain 
sharp  corners  and  steps.  These  sharp  geometric  variations  make  it  extremely 
difficult  to  generate  body-conforming  grids  while  preserving  the  sharp  cor¬ 
ners.  The  grid  lines  are  wrapped  around  the  corners  and,  in  many  cases,  such 
wrap  around  grids  are  skewed  near  these  corners  and  steps.  Using  such  grids 
introduces  geometric  errors  and  sometimes  leads  to  loss  in  both  the  computa¬ 
tional  efficiency  and  accuracy.  In  this  report  we  develop  and  apply  a  flow- 
field  blanking  procedure  which  allows  computation  of  practical  flows  of  inter¬ 
est  with  no  geometric  error  since  it  models  the  corners  and  steps  exactly. 

To  avoid  geometric  errors  one  can  blank  out  the  flowfield  in  specific 
regions  in  the  computational  domain.  Examples  where  such  blanking  can  be 
useful  are  shown  In  Figure  1.  Continuous  straight  line  grids  can  be  used  for 
these  cases  and  the  hatched  regions  are  the  ones  where  the  flowfield  is  to  be 
blanked  out.  This  procedure,  thus,  preserves  the  sharp  corners  and  steps.  In 
addition  to  zeroing  out  the  flowfield  Inside  the  hatched  regions,  additional 
changes  must  be  made  in  the  boundary  conditions  and  the  computational  algo¬ 
rithm  near  these  surfaces.  These  changes  are  described  in  a  later  section. 
This  technique  can  be  tested  with  the  simple  problan  of  flow  over  a  rotating 
band.  The  rotating-band,  which  is  a  protuberance  on  the  artillery  shell, 
imparts  spin  to  a  shell  during  launch.  However,  it  does  contribute  a  small 
unwanted  drag  in  free  flight.  A  schematic  of  the  rotating-band  flowfield  is 
shown  in  Figure  2.  It  shows  the  expected  recirculation  regions  in  front  of 
and  behind  the  band  and  the  associated  compressions  and  expansion  waves.  A 
numerical  solution  is  obtained  for  this  problem  at  M  =  3.0  and  a  =  0. 


COMPUTATIONAL  TECHNIQUE 


GOVERNING  EQUATIONS 


The  complete  set  of  time-dependent  generalized  axi symmetric  thin- layer 
Navier-Stokes  equations  is  solved  numerically  to  obtain  a  solution  to  this 
problem.  The  numerical  technique  used  is  an  implicit  finite-difference 
scheme.  Although  time-dependent  calculations  are  made,  the  transient  flow  is 
not  of  primary  interest  at  the  present  time.  The  steady  flow  is  the  desired 
result  which  is  obtained  in  a  time  asymptotic  fashion. 

The  azimuthal -i nvari ant  (or  generalized  axi symmetri c)  tnin-layer  Navier- 
Stokes  equations  for  curvilinear  coordinates  £,  n  and  t  can  be  written  as:1 


C(x,y,z,t)  is  the  longitudinal  coordinate 
n(y,z,t)  is  the  circumferential  coordinate 
c(x,y,z,t)  is  the  near  normal  coordinate 
t  is  the  time 
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The  velocities 


U  =  5t  +  Cxu  +  £yv  *  Czw 

V  =  nt  *  nxu  +  nyv  +  n2w  (2) 

W  =  Ct  +  CXU  +  cyv  +  ;2w 

represent  the  contravariant  velocity  components. 

The  Cartesian  velocity  components  (u,  v,  w)  are  nondimensionalized  with 
respect  to  a^  (free  stream  speed  of  sound).  The  density  (p)  is  referenced 

To  Pa,  and  total  energy  (e)  to  P.a*2.  The  local  pressure  is  determined  using 

the  equation  of  state, 

p  *  (y  -  l)[e  -  0.5p(u2  +  v2  +  w2)]  (3) 


where  y  is  the  ratio  of  specific  heats. 

While  Equation  (1)  contains  only  two  spatial  derivatives,  it  retains  all 
three  momentum  equations,  thus  allowing  a  degree  of  generality  over  the  stan¬ 
dard  axisymmetric  equations.  In  particular,  the  circumferential  velocity  is 
not  assumed  to  be  zero,  thus  allowing  computations  for  spinning  projectiles  or 
swirl  flow  to  be  accomplished. 

2.  COMPUTATIONAL  ALGORITHM 

The  azimuthal-invariant,  thin-layer  Navi er-Stokes  equations  are  solved 
using  an  implicit  approximate  factorization  finite  difference  scheme  in  delta 
form.6  An  implicit  method  was  chosen  because,  for  viscous  flow  problems,  it 
permits  a  time  step  much  greater  than  that  allowed  by  explicit  schemes.  The 
Beam-Wanning  implicit  algorithm  has  been  used  in  various  applications1”9  for 
the  equations  in  general  curvilinear  coordinates.  The  algorithm  is  first-order 
accurate  in  time  and  second-  or  fourth-order  accurate  in  space.  The  equations 
are  factored  (spatially  split),  which  reduces  the  solution  process  to  one¬ 
dimensional  problems  at  a  given  time  level.  Central  difference  operators  are 
employed  and  the  algorithm  produces  block  tridiagonal  systems  for  each  space 
coordinate.  The  main  computational  work  is  contained  in  the  solution  of  these 
block  tridiagonal  systems  of  equations.  For  the  computation  of  turbulent 
flows,  the  two-layer  algebraic  Bal dwi n-Lomax  turbulence  model10  is  used. 

3.  FINITE-DIFFERENCE  EQUATIONS 

The  implicit,  approx imately  factored  algorithm  developed  by  Bean- 
Warming6  has  the  form: 
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(4) 


[I  +  h<5  An  +  D  ]  [I  +  h<5  Cn  -  hRe'1  6  J'1  MnJ  +  0  ]Aqn 

S  S  C  £  S 


=  -h[6^En  +  6?Gn  -  Re'1  5^Sn  +  Hn  +  ] 

where  the  explicit  fourth-order  dissipation  is: 

(<♦) 

0  =  -eeAtJ_1[(75Ac)2  +  (7^?)2]J  qn 

and  the  implicit  second-order  dissipation  terms  are: 

( 2 ) 

=  -  Ei AtJ  1(7?Ac)J 

(2) 

D;  =  -  £iAtJ  1(7cAc)J  . 


The  fourth-order  explicit  dissipation  is  used  to  control  non-linear  instabil¬ 
ities  whereas  the  implicit  dissipation  is  included  to  stabilize  the  explicitly 
treated  fourth-difference  terms.  The  parameter  eg  is  0(1)  and  the  parameter 

*  3  H  A  ^  q 

e,  is  two  to  three  times  e  .  The  Jacobian  matrices  A  =  — *  ,  C  *  — along  with 

.  A  3q  3q 

coefficient  matrix  M  obtained  from  linearization  of  S  are  described  in  detail 
in  Reference  8. 

To  suppress  high  frequency  components  that  appear  in  regions  containing 
severe  pressure  gradients,  e.g.,  shocks  or  stagnation  points,  a  switching 
dissipation  model  is  used.  This  switching  model  is  similiar  to  the  model  used 
by  Pulliam9  and  uses  a  fourth-order  dissipation  in  smooth  regions  and  switches 
to  a  second-order  dissipation  in  regions  containing  high  pressure  or  density 

U) 

gradients.  The  dissipation  term  n  on  the  right  hand  side  of  Equation  (4) 
can  he  written  in  this  model  as: 


— ■  I  I  I  I  l-^yl  5  J  q  -  6  eg  5A7  J  q] 


(S) 


where  the  first  term  is  the  second-order  dissipation  and  the  second  term  con¬ 
tains  the  fourth-order  dissipation.  The  coefficients  ed  and  eg  are  the  asso¬ 
ciated  coefficients  for  the  second-order  and  fourth-order  dissipation,  respec¬ 
tively.  The  coefficient  ed  is  fifty  to  hundred  times  and  A  and  7  are  the 

one-sided  forward  and  backward  finite-difference  operators.  Mote  tnat  the 
fourth-order  dissipation  is  non-linear  in  that  the  coefficient  is  not  a  con¬ 
stant  and  is  scaled  hy  spectral  radius  IjA^JI.  The  two  terms  in  Equation  [S' 
are  of  the  form  6a6e  where: 


(6  a  5  B)j  =  (- 


°j  +  aj-l 


- - I)  (3-  .  -  B-)  +  (— - 

2  1  VPj+i  V  v  2 


)  (Bj  -  Bj.i) 


Fourth-order  dissipation  is  used  if  eg  >  and  the  dissipation  is 

switched  to  second-order  if  eg  <  ^7^  *  The  Pressure  gradient  is  used  in 

the  normal  direction  in  this  switching  control  whereas  density,  as  shown  in 

Equation  (5),  is  used  in  the  longitudinal  direction.  In  addition,  a  space 
varying9  At  procedure  is  used  where  the  time  step  used  is  given  as: 


1  ♦  ZT 

where  J  is  the  Jacobian  of  the  transformation  and  (At)  ^  is  a  reference  time 
step. 

4.  FLOWFIELf)  BLANKING 

The  idea  is  to  avo1d  geometric  errors  that  may  arise  from  wrap  around 
grids.  Instead,  we  use  straight  line  grids  as  shown  schematically  in  Figure 
3.  For  the  rotating  band  problem,  the  zone  ABCD  is  part  of  the  body  and  the 
flowfield  in  this  zone  must  be  blanked  out  in  the  computational  domain.  As 
shown  in  Figure  3,  the  sharp  corners  and  steps  ahead  of  and  behind  the  band 
are  preserved  and  no  approximation  is  made.  It  is  also  necessary  to  apply 
boundary  conditions  on  the  zonal  surfaces  AB,  BC  and  CD.  The  no-slip  boundary 
conditions  are  used  at  these  boundaries  along  with  zero  gradients  for  pressure 
and  density.  In  addition,  at  neighboring  points  to  these  boundaries,  we  use 
second-order  spatial  difference  and  smoothing.  The  block  tri diagonal  matrix 
structure  has  been  modified  for  continuous  integration  sweeps  through  such 
zones.  For  example,  the  block  tridiagonal  matrix  in  the  £  direction  takes  the 
following  form  (after  setting  =  0  to  simplify  the  illustration) 


"AJMAX-2  1 


AqJMAX-l 


*HSJMAX-1 
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Here  A's  denote  the  quantity  A  and  I  is  a  5  x  5  identity  matrix.  Note  the 

appearance  of  the  uncoupled  block  tridiagonals  between  J  =  J1  and  J2  corre¬ 
sponding  to  lines  A8  and  DC,  respecti vely.  The  rows  at  J1  and  J2  are  particu¬ 
larly  simple  because  boundary  conditions  are  updated  explicitly  at  the  end  of 
inversions.  All  the  changes  described  in  this  section  were  easily  implemented 
in  a  modular  fashion  into  an  existing  code  for  projectile  flow  computations. 
One  simply  fills  the  block  tridiagonal  matrix  ignoring  the  zone.  Elements  in 
the  rows  inside  the  zone  are  then  overloaded  as  shown  above.  The  flowfield 
blanking  affects  the  block  tridiagonal  matrix  in  the  c  direction  similarly. 
Although,  we  have  only  one  zone  for  the  rotating  band  case,  changes  have  been 
made  in  the  code  to  blank  out  multiple  zones. 


III.  RESULTS 

All  the  numerical  computations  were  made  at  =  3.0  and  a  =  0.  The 

projectile  configuration  with  the  rotating  band  which  was  used  in  this  study 
is  shown  in  Figure  4.  This  model  is  a  cone-cyl inder  configuration  with  a 
13.1°  cone  angle.  The  band  height  is  .04  0  and  the  width  is  .505  D.  The  same 
model  was  used  in  the  experiments1 1  which  were  conducted  in  the  US  Army  Chem¬ 
ical  Research  Development  and  Engineering  Center's  Supersonic  Wind  Tunnel. 
Surface  pressure  measurements  have  been  made  ahead  of  and  behind  the  band 
which  are  used  to  compare  with  the  numerical  results. 

Since  the  freestream  flow  is  supersonic,  the  space  marching  Parabolized 
Navier-Stokes  code3  was  used  to  compute  the  solution  over  the  forebody  of  the 
projectile  (See  Figure  5).  This  generated  a  solution  at  a  station  30  band 
heights  ahead  of  the  band  which  was  then  used  as  an  upstream  boundary  con¬ 
dition  for  the  computation  of  the  flowfield  containing  the  rotating  band.  For 
this  part  of  the  flowfield  which  includes  the  band,  the  unsteady  or  time- 
dependent  Navier-Stokes  computational  technique  described  earlier  was  used. 
Such  composite  solution  technique  allowed  a  large  number  of  grid  points  to  be 
used  in  the  vicinity  of  the  band. 

The  computati onal  grid  used  for  the  numerical  calculations  is  shown  in 
Figure  5.  It  consists  of  130  points  in  the  longitudinal  direction  and  61 
points  in  the  normal  direction.  The  grid  points  are  clustered  near  the  sur¬ 
face  of  the  cylindrical  part  with  a  minimum  spacing  of  .00002  0.  The  resolu¬ 
tion  of  grid  points  on  the  top  of  the  band  is  not  as  fine.  Grid  points  in  the 
longitudinal  direction  are  clustered  near  the  upstream  and  downstream  corners 
of  the  rotating  band  where  appreciable  changes  in  the  flow  variables  are  ex¬ 
pected.  In  Figure  6,  the  grid  lines  inside  the  band  are  omitted  to  show  the 
position  of  the  band;  however,  in  the  actual  grid  used  in  the  computations, 
there  are  continuous  grid  lines  inside  the  band  and  those  are  the  lines  where 
the  flowfield  blanking  procedure  is  used. 

For  comparison  purposes,  a  numerical  solution  is  first  obtained  for  flow 
over  the  cylindrical  part  of  the  projectile  without  the  rotating  band  at  = 

3.1  and  u  =  1.  The  computed  surface  pressure  coefficient  is  plotted  in  Figure 
7  as  a  function  of  longitudinal  position.  The  computed  result  is  in  very  good 
agreement  with  experimental  data.11 
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Numerical  results  obtained  for  the  rotating  band  case  are  presented 
next.  Figure  3a  shows  the  velocity  vector  field  in  front  of  the  band  and  as 
expected,  it  shows  the  reci rcul atory  flow  in  that  region.  As  shown  in  this 
figure,  the  flow  seems  to  accelerate  as  the  corner  of  the  band  is  approached. 
Figure  3b  shows  the  velocity  vectors  behind  the  band.  The  flow  expands  at  the 
corner  of  the  band.  A  reci rcul ation  region  can  be  observed  clearly.  Figures 
9a  and  9b  show  the  stream  function  contours  ahead  of  and  behind  the  band,  re¬ 
spectively.  The  recirculatory  flow  regions  can  be  clearly  seen  in  these  fig¬ 
ures.  The  reverse  flow  region  extends  about  four  band  heights  ahead  of  the 
band  and  the  reattachment  point  is  less  than  a  quarter  of  the  height  of  the 
band  from  the  corner.  The  size  of  the  recirculation  bubble  behind  the  band  is 
a  little  smaller  than  the  one  ahead  of  the  band.  The  flow  seems  to  separate 
slightly  below  the  band  corners  and  reattaches  about  3.5  band  heights  down¬ 
stream.  Figure  10  shows  the  pressure  contours  for  this  case.  One  can  also 
see  a  separation  shock  wave  ahead  of  the  band.  The  shock  wave  is  located  just 
ahead  of  the  flow  separation  region.  The  strong  flow  expansions  at  both  the 
band  corners  can  be  clearly  seen.  The  expansions  at  the  downstream  corner  is 
followed  by  a  recompression  shock.  The  surface  pressure  coefficient  for  the 
band  case  is  shown  in  Figure  11  as  a  function  of  the  axial  position.  The 
solid  line  is  the  computed  result,  the  dashed  line  is  the  result  obtained  for 
the  case  without  the  band  and  the  circles  are  the  experimental  data  for  the 
band.  There  is  a  considerable  change  in  the  pressure  due  to  the  presence  of 
the  band.  The  sharp  rise  in  pressure  ahead  of  the  band  is  associated  with  the 
shock  wave  which  actually  precedes  the  separation  point  of  the  boundary  layer 
flow.  The  flow  then  expands  at  the  corner  and  pressure  drops.  No  significant 
change  in  pressure  occurs  on  the  top  of  the  band.  At  the  backward  step  of  the 
band,  the  flow  expands  again  which  results  in  the  sharp  decrease  in  the  pres¬ 
sure.  This  is  followed  by  a  more  gradual  return  to  the  ambient  pressure  down¬ 
stream.  The  computed  surface  pressure  is  in  good  agreement  with  the  experi¬ 
mental  data  measured  ahead  of  and  behind  the  band.  The  small  discrepancy 
found  in  the  comparison  could  be  due  to  the  turbulence  model  used. 


IV.  CONCLUDING  REMARKS 

The  flavier-Stokes  computational  technique  has  been  used  in  conjunction 
with  a  flowfield  blanking  procedure  for  numerical  simulation  where  the  sharp 
corners  and  steps  exactly  modeled,  thereby,  avoiding  any  possible  source  of 
geometric  errors.  This  procedure  has  been  applied  to  the  flow  over  a  rotating 
band  at  supersonic  speed. 

Computed  results  have  been  obtained  for  =  3.0  and  o  =  0  and  compared 

with  available  experimental  data.  The  results  show  the  reci  rcul  ati  on  region 
both  ahead  of  and  behind  the  rotating-band  as  well  as  the  associated 
compression  and  expansion  waves.  The  computed  surface  pressures  for  both 
cases,  with  and  without  the  band,  are  in  fairly  good  agreement  with  experi¬ 
mental  data.  The  present  numerical  procedure  is  simple  to  use  and  seems  to 
predict  the  flowfield  correctly.  Further  work  is  needed  to  extend  this 
technique  to  predict  three  dimensional  flow  fields.  In  addition,  a  parametric 
study  will  be  conducted  in  future  to  predict  the  effect  of  the  rotati ng-band 
on  the  aerodynamic  coefficients  for  artillery  shell. 
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LIST  OF  SYMBOLS 


a  =  speed  of  sound 

Cp  =  specific  heat  at  constant  pressure 

Cp  =  pressure  coefficient,  2(p0Oa2p  -  pj/p.u* 

D  =  body  diameter 

e  =  total  energy  per  unit  volume/p^a* 

E,  F,  q  =  flux  vector  of  tranformed  Navi er- Stokes  equations 
H  =  n-invariant  source  vector 

J  =  Jacobian  of  transformati on 

M  =  Mach  number 

p  =  pressure/ p^a^ 

Pr  =  Prandtl  number,  p  e  /< 

ao  Q  ao 

R  =  body  radius 

Re  =  Reynolds  number,  P^ajl/p^ 

S  =  viscous  flux  vector 

t  =  physical  time 

u,v,w  =  Cartesian  velocity  component s/a^ 

U  , V , W  =  contravari ant  velocity  component s/a^ 

x,y,z  =  physical  Cartesian  coordinates 

a  =  angle  of  attack 

y  =  ratio  of  specific  heats 

<  =  coefficient  of  thermal  conductivity/^ 

p  =  coefficient  of  vi scosi ty/p^ 

L,n,c  -  transformed  coordinates  in  axial,  circumferential  and  radial 
di reef  1 ons 

U  =  density/p^ 

forward  difference 
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LIST  OF  SYMBOLS  (Continued) 


backward  difference 

central  difference 

transformed  time 

implicit  smoothing  coefficient 

second  order  dissipation  coefficient 

fourth  order  dissipation  coefficient 

critical  value 

longitudinal  direction 
identity  matrix 

free  stream  conditions  for  corresponding  dimensional  quantity 
streamwise  direction 


normal  direction 
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